$ontext

   An example of Benders Decomposition on fixed charge transportation
   problem bk4x3. This formulation has a few refinements to speed
   up convergence.

   Optimal objective in reference : 350.

   Erwin Kalvelagen, December 2002

   See:
   http://www.in.tu-clausthal.de/~gottlieb/benchmarks/fctp/
   http://www.amsterdamoptimization.com/pdf/benders.pdf


$offtext


set i 'sources' /i1*i4/;
set j 'demands' /j1*j3/;

parameter supply(i) /
   i1 10
   i2 30
   i3 40
   i4 20
/;


parameter demand(j) /
   j1 20
   j2 50
   j3 30
/;


table c(i,j) 'variable cost'
       j1    j2    j3
i1    2.0   3.0   4.0
i2    3.0   2.0   1.0
i3    1.0   4.0   3.0
i4    4.0   5.0   2.0
;



table f(i,j) 'fixed cost'
        j1     j2     j3
i1    10.0   30.0   20.0
i2    10.0   30.0   20.0
i3    10.0   30.0   20.0
i4    10.0   30.0   20.0
;

*
* check supply-demand balance
*
scalar totdemand, totsupply;
totdemand = sum(j, demand(j));
totsupply = sum(i, supply(i));
abort$(abs(totdemand-totsupply)>0.001) "Supply does not equal demand.";

*
* for big-M formulation we need tightest possible upperbounds on x
*
parameter xup(i,j) 'tight upperbounds for x(i,j)';
xup(i,j) = min(supply(i),demand(j));


*--------------------------------------------------------------------
* standard MIP problem formulation
*--------------------------------------------------------------------

variables
   cost   'objective variable'
   x(i,j) 'shipments'
   y(i,j) 'on-off indicator for link'
;
positive variable x;
binary variable y;

equations
   obj      'objective'
   cap(i)   'capacity constraint'
   dem(j)   'demand equation'
   xy(i,j)  'y=0 => x=0'
;

obj..        cost =e= sum((i,j), f(i,j)*y(i,j) + c(i,j)*x(i,j));
cap(i)..     sum(j, x(i,j)) =l= supply(i);
dem(j)..     sum(i, x(i,j)) =g= demand(j);
xy(i,j)..    x(i,j) =l= xup(i,j)*y(i,j);

model fscp /obj,cap,dem,xy/;

option optcr=0;
option limrow=0;
option limcol=0;


*---------------------------------------------------------------------
* Benders Subproblem
*---------------------------------------------------------------------

variable z 'objective variable';

positive variables
   u(i) 'duals for capacity constraint'
   v(j) 'duals for demand constraint'
   w(i,j) 'duals for xy constraint'
;

equations
   subobj          'objective'
   subconstr(i,j)  'dual constraint'
;

* to detect unbounded subproblem
scalar unbounded /1.0e6/;
z.up = unbounded;

subobj..  z =e= sum(i, -supply(i)*u(i)) + sum(j, demand(j)*v(j))
                + sum((i,j), -xup(i,j)*y.l(i,j)*w(i,j))
                 ;

subconstr(i,j)..  -u(i) + v(j) - w(i,j) =l= c(i,j);

model subproblem /subobj, subconstr/;
* reduce output to listing file:
subproblem.solprint=2;
* speed up by keeping GAMS in memory:
subproblem.solvelink=2;

*---------------------------------------------------------------------
* Benders Modified Subproblem to find unbounded ray
*---------------------------------------------------------------------

variable dummy 'dummy objective variable';

equations
   modifiedsubobj          'objective'
   modifiedsubconstr(i,j)  'dual constraint'
   edummy;
;

modifiedsubobj..
    sum(i, -supply(i)*u(i)) + sum(j, demand(j)*v(j))
           + sum((i,j), -xup(i,j)*y.l(i,j)*w(i,j)) =e= 1;

modifiedsubconstr(i,j)..
    -u(i) + v(j) - w(i,j) =l= 0;

edummy.. dummy =e= 0;

model modifiedsubproblem /modifiedsubobj, modifiedsubconstr, edummy/;
* reduce output to listing file:
modifiedsubproblem.solprint=2;
* speed up by keeping GAMS in memory:
modifiedsubproblem.solvelink=2;


*---------------------------------------------------------------------
* Benders Relaxed Master Problem
*---------------------------------------------------------------------

set iter /iter1*iter50/;

set cutset(iter) 'dynamic set';
cutset(iter)=no;
set unbcutset(iter) 'dynamic set';
unbcutset(iter)=no;


variable z0 'relaxed master objective variable';
equations
   cut(iter)           'Benders cut for optimal subproblem'
   unboundedcut(iter)  'Benders cut for unbounded subproblem'
   ycon1(j)            'extra condition on y'
   ycon2(i)            'extra condition on y'
;

parameters
   cutconst(iter)     'constant term in cuts'
   cutcoeff(iter,i,j)
;

cut(cutset).. z0 =g= sum((i,j), f(i,j)*y(i,j))
                      + cutconst(cutset)
                      + sum((i,j), cutcoeff(cutset,i,j)*y(i,j));
unboundedcut(unbcutset)..
                cutconst(unbcutset)
                + sum((i,j), cutcoeff(unbcutset,i,j)*y(i,j)) =l= 0;

ycon1(j)..        sum(i, y(i,j)*supply(i)) =g= demand(j);
ycon2(i)..        sum(j, y(i,j)*demand(j)) =g= supply(i);


model master /cut,unboundedcut,ycon1,ycon2/;
* reduce output to listing file:
master.solprint=2;
* speed up by keeping GAMS in memory:
master.solvelink=2;
* solve to optimality
master.optcr=0;


*---------------------------------------------------------------------
* Benders Decomposition Initialization
*---------------------------------------------------------------------


display "--------------------- BENDERS ALGORITHM ----------------------------";

scalar UB 'upperbound' /INF/;
scalar LB 'lowerbound' /-INF/;

model feasy /edummy,ycon1,ycon2/;
* reduce output to listing file:
feasy.solprint=2;
* speed up by keeping GAMS in memory:
feasy.solvelink=2;
solve feasy minimizing dummy using mip;

display "Initial values",y.l;


*---------------------------------------------------------------------
* Benders Algorithm
*---------------------------------------------------------------------

scalar converged /0/;
scalar iteration;
scalar bound;
parameter ybest(i,j);
parameter log(iter,*) 'logging info';


loop(iter$(not converged),

*
* solve Benders subproblem
*
   solve subproblem maximizing z using lp;

*
* check results.
*

   abort$(subproblem.modelstat>=2) "Subproblem not solved to optimality";

*
* was subproblem unbounded?
*

   if (z.l+1 < unbounded,

*
* no, so update upperbound
*

      bound = sum((i,j), f(i,j)*y.l(i,j)) + z.l;
      if (bound < UB,
          UB = bound;
          ybest(i,j) = y.l(i,j);
          display ybest;
      );


*
* and add Benders' cut to Relaxed Master
*
      cutset(iter) = yes;

   else

*
* solve modified subproblem
*

     solve modifiedsubproblem maximizing dummy using lp;

*
* check results.
*

     abort$(modifiedsubproblem.modelstat>=2)
            "Modified subproblem not solved to optimality";


*
* and add Benders' cut to Relaxed Master
*
      unbcutset(iter) = yes;
   );


*
* cut data
*
   cutconst(iter) = sum(i, -supply(i)*u.l(i)) + sum(j, demand(j)*v.l(j));
   cutcoeff(iter,i,j) = -xup(i,j)*w.l(i,j);

*
* solve Relaxed Master Problem
*

   option optcr=0;
   solve master minimizing z0 using mip;

*
* check results.
*

   abort$(master.modelstat=4) "Relaxed Master is infeasible";
   abort$(master.modelstat>=2) "Masterproblem not solved to optimality";

*
* update lowerbound
*

   LB = z0.l;

   log(iter,'LB') = LB;
   log(iter,'UB') = UB;

   iteration = ord(iter);
   display iteration,LB,UB;

   converged$( (UB-LB) < 0.1 ) = 1;
   display$converged "Converged";

);

display log;

abort$(not converged) "No convergence";

*
* recover solution
*
y.fx(i,j) = ybest(i,j);
fscp.solvelink=2;
fscp.solprint=2;
solve fscp minimizing cost using rmip;
abort$(fscp.modelstat<>1) "final lp not solved to optimality";

display "Benders solution",y.l,x.l,cost.l;